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We analyze the spatial correlation structure of the spin density of an electron gas in the vicinity 
of an antiferromagnetically-coupled Kondo impurity. Our analysis extends to the regime of spin- 
anisotropic couplings, where there are no quantitative results for spatial correlations in the literature. 

We use an original and numerically exact method, based on a systematic coherent-state expansion 
of the ground state of the underlying spin-boson Hamiltonian. It has not yet been applied to the 
computation of observables that are specific to the fermionic Kondo model. We also present an 
important technical improvement to the method, that obviates the need to discretize modes of the 
Fermi sea, and allows one to tackle the problem in the thermodynamic limit. As a result, one can 
obtain excellent spatial resolution over arbitrary length scales, for a relatively low computational 
cost, a feature that gives the method an advantage over popular techniques such as the Numerical 
and Density-Matrix Renormalization Groups. We find that the anisotropic Kondo model shows rich 
universal scaling behavior in the spatial structure of the entanglement cloud. First, SU(2) spin- 
symmetry is dynamically restored in a finite domain in parameter space in vicinity of the isotropic 
line, as expected from poor man’s scaling. More surprisingly, we are able to obtain in closed 
analytical form a set of different, yet universal, scaling curves for strong exchange asymmetry, 
which are parametrized by the longitudinal exchange coupling. Deep inside the cloud, i.e. for 
distances smaller than the Kondo length, the correlation between the electron spin density and the 
impurity spin oscillates between ferromagnetic and antiferromagnetic values at the scale of the Fermi 
wavelength, an effect that is drastically enhanced at strongly anisotropic couplings. Our results also 
provide further numerical checks and alternative analytical approximations for the Kondo overlaps 
that were recently computed by Lukyanov, Saleur, Jacobsen, and Vasseur [Phys. Rev. Lett. 114, 

080601 (2015)] . 

PACS numbers: 73.40.Gk, 72.10.Fk 


I. INTRODUCTION 

The spin 1/2 Kondo model describes a localized mag¬ 
netic moment interacting with an electron gas, via an 
antiferromagnetic exchange coupling .^1 Despite a long his¬ 
tory, and even an exact solution, it has not yet sur¬ 
rendered all its secrets. It is well established that the 
ground state is a spin singlet in which the impurity spin 
is quenched by the electron gas, and the spatial region 
where the electron gas is correlated with the impurity is 
referred to as the Kondo screening cloud.^^ Even for an 
isotropic system, its precise spatia^rohle is not known 
analytically, except asymptotically, 121 and it is only in the 
past few years that it has been calculated numerically.^^®] 
Despite wide-ranging proposals, eluded direct 
measurement, partly because of the difficulty in measur¬ 
ing spin correlations. 

In the isotropic case, the screening cloud is charac¬ 
terized by the ground state correlation function X(x) = 

4 • 5®*(a;)^ where is the impurity spin opera¬ 

tor, and 5®'(a:) is the electron spin density at x. When 
the Kondo temperature is much lower than the Fermi 
energy, the screening cloud can be decomposed into a 


forward scattering contribution Xq{x) and a backscat- 
tering contribution X 2 kp{x), so that X{x) = Xq{x) -|- 
cos{2kpx)X2kF{x), where the two functions Xq{x) and 
X 2 kp (x) vary slowly on the scale of the Fermi wave¬ 
length 2TT/kp. In the scaling regime and for spin-isotropic 
Kondo exchange, the profil e of t he screening cloud dis¬ 
plays a universal line shape,^i2E^ which is dependent on 
the value of Kondo coupling J only through an emer¬ 
gent length f that is inversely proportional to the Kondo 
temperaturel ^^ l ^^ l To be specihc, if Xk{x) and X'^{x), 
with k S {0,2^^}, are correlation functions correspond¬ 
ing to different values J and J' of the Kondo coupling, 
and ^ and are the corresponding Kondo lengths, then 

X/(x) = ^Xk{f,'x/£,), (1) 

for all X » 2'Klkp. 

Recently, we have realized that it may be possible 
to measure the longitudinal forward scattering [k = 0) 
component of the screening cloud (a precise definition 
is given below) in a chain of tunnel-coupled supercon¬ 
ducting islands.^ It turns out that the Hamiltonian de¬ 
scribing the charge and phase degrees of freedom in this 
system is equivalent to the spin-anisotropic Kondo model 
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in one dimension,^ described by two coupling constants: 
The ^-components of the impurity and electron spins cou¬ 
ple with a strength Jy, while the components perpendic¬ 
ular to the z-axis couple with a different strength J_l . A 
detailed study of the screening cloud in the anisotropic 
Kondo model has to our knowledge not been performed 
before, and is therefore timely. 

In the anisotropic case, there are four correlation func¬ 
tions of interest, namely AjJ(a:), that measures the cor¬ 
relation between z-components of the impurity and elec¬ 
tron spins (both in the forward k = Okp and backward 
k = 2kF scattering channels), and X^{x), that measures 
correlations between components perpendicular to the z 
axis (in each channel). Clearly, an anisotropic Kondo in¬ 
teraction, characterized by two coupling constants Jy and 
J_L, will affect the universal scaling picture non-trivially. 

The universal scaling of the isotropic model has been 
confirmed using the Numerical Renormalization Gr oup 
(NRG), a popular method to study the Kondo model.l^ 
However, there is still room to improve the accuracy of 
existing results. For instance, it is predicted analyticalljPl 
that backscattering dominates forward scattering inside 
the screening cloud. This produces oscillations from fer¬ 
romagnetic to antiferromagnetic correlations between the 
impurity and the electron gas, on the scale of half the 
Fermi wavelength. The first numerical calculation of the 
universal scaling functions were reported in Ref. [S) How¬ 
ever, the predominance of the backscattering component 
over the forward scattering component inside the cloud 
was not resolved. This was probably due to the fact that 
results were obtained at relatively large Kondo temper¬ 
atures, so that there was not a sufficient separation of 
scales between the size ^ of the cloud, and the short dis¬ 
tance ultraviolet cut-off scale. In a more recent work,l^ 
the NRG method of Ref. was refined, and the alter¬ 
nation of ferro- and antiferromagnetic correlations inside 
the cloud clearly be seen. However, results were only pre¬ 
sented for distances up to 10 Fermi wavelengths from the 
impurity. At these scales, the correlation functions show 
non-universal modulations. 

In this Article, we perform numerical calculations that 
are sufficiently accurate to investigate the universal scal¬ 
ing behavior in regimes ranging from isotropic to strongly 
anisotropic Kondo couplings. At the same time, we 
can clearly resolve the predominace of the backscat¬ 
tering component over the forward scattering compo¬ 
nent inside the cloud. In precise terms, we investigate 
the the following questions: Let Xl{x) and {Xiy{x), 
with j e {||,T} and k G {0,2kF}, be correlation func¬ 
tions at distinct values (JL,Jy) and (J(|_,J|P. Under 
which conditions are there constants Ai and A 2 such that 
inYi x) = AiX^(A 2 x), for all x ^ 2FlkF-, and what are 
the line shapes of these universal scaling curves? Our 
general numerical findings show that scaling is typically 
obeyed for fixed values of Jy. In other words, correla¬ 
tors at the same Jy but different Jx can be scaled onto 
each other. This statement holds as long as Jx ^ I 
(in dimensionless units of the inverse density of states). 


while for transverse couplings of order one, the Kondo 
temperature becomes comparable to the Fermi energy, 
and universality is lost. For Jy ^ Jx, we find that the 
line shapes of the cloud correlation functions acquire a 
Jy -dependence that cannot be scaled away. However, as 
Jy approaches Jx, the Jy dependence rapidly becomes 
very weak. One thus find a sizable region of parame¬ 
ter space, located around the strict spin-isotropic line 
J\\ = J± ^ 1, in which the universal curves have lit¬ 
tle discernible Jy dependence. In this region, correlators 
calculated at Jx Jx and Jy ^ Jy can to a very good 
approximation be scaled onto the unique universal curve 
of the isotropic model. These conclusions are consistent 
with poor man’s scaling argument^^ at small Jx, as we 
will explain in detail. 

Our work also elaborates on a different but related 
topic, that was raised very recently by Lukyanov et al. 
in Ref. [5TJ where the overlap between two Kondo wave 
functions with different Kondo couplings was computed 
analytically using integrability techniques, and numeri¬ 
cally with the Density Matrix Renormalization Group. 
We complement these interesting results, by providing 
a simpler (approximate) analytical expression for the 
Kondo overlaps in the spin-anisotropic limit. We also 
show that our numerical method reproduces the full an¬ 
alytical result, as it should. 

Regarding methodology, our approach combines both 
analytics and numerically exact calculations. In the case 
of strong spin-anisotropy, Jy ^ Jx, we provide asymp¬ 
totically exact formulas, that were not discussed in pre¬ 
vious literature, for the fermionic Kondo cloud. For a 
general choice of parameters, we employ a recent ap¬ 
proach using a coherent state expansion of the wave 
function,l22E31 that has not yet been applied to the com¬ 
putation of purely fermionic observables that pertain to 
the Kondo model. This method relies on a variational 
Ansatz, formulated in the language of the spin-boson 
model, which describes a two-lev el sy stem coupled to 
a multi-mode ohmic bosonic bath.l^^l^SI This spin-boson 
model is known to be equivalent to the anisotropic Kondo 
model,^ with b osoniz ation providing an exact mapping 
between the twol ^^ l ^^ l We exploit this mapping to apply 
the coherent-state expansion to the Kondo model. In ad¬ 
dition, the flexible structure of the Ansatz allows one to 
add progressively more contributions to it, so that the 
method converges rapidly to the exact ground state with 
arbitrary accuracy. 

Historically, the syst emati c coherent state expansion 
proposed by Bera et builds on semin al wo rks of 

Emery and Luther,!^ and of Silbey and Harris ,^221211 in the 
language of the spi n-bos on model, and on independent 
works by Anderson ^nd Bergmann and Zhan^^ in 
the language of the Kondo model. The approach of 
Emery and Luther, Silbey and Harris, and Anderson cor¬ 
responds to the lowest order approximation, involving a 
single coherent state, and provides quantitatively accu¬ 
rate results only in the limit Jy 3> Jx- Bergmann and 
Zhang took a step in the direction of a two coherent- 
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state Ansatz by adding an extra term to the lowest order 
approximation, but constrained its form sub-optimally 
compared to the full solution with two coherent states. 
We stress that this type of approach differs from the 
well-known method pioneered by Yosida,^ that takes as 
starting point a state in which a single electron binds 
into a singlet with the impurity, while the rest of the 
Fermi sea is unaffected. The wave function of the bound 
electron is chosen to maximize the binding energy. The 
relationship between the Yosida approximation and the 
single-coherent state Ansatz is reminiscent of the rela¬ 
tionship in superconductivity between the single Cooper 
pair on top of a Fermi sea on the one hand, and the full 
BCS wave function on the other. Yosida’s approach gives 
qualitatively correct results, but even when the effect of 
additional particle-hole excitations are included it 
cannot yield arbitrarily accurate results for the correla¬ 
tions inside the Kondo screening cloud. 

In this Article, we also contribute an important in¬ 
novation to the coherent-state expansion methodology. 
In previous implementations, the bosonic bath was lim¬ 
ited to a finite number of modes, because the number of 
variational parameters in the Ansatz scaled linearly with 
the number of bath modes. Here, we first partially solve 
the variational problem analytically, so that the number 
of remaining variational parameters that have to be de¬ 
termined numerically is independent of the number of 
bath oscillators. This allows us to take the thermody¬ 
namic limit before we optimize the energy numerically, 
and probe arbitrarily large distances. In the Kondo lan¬ 
guage, this means that we work directly with an infinite 
conductor. 


With NRG, probing large length scales comes at the 
cost of losing information about shorter lengths scales. 
In order to maintain good spatial resolution, one has to 
introduce a fictitious impurity at the position where cor¬ 
relators are evaluated. Each position considered then 
requires the solution of a given two-channel Kondo im¬ 
purity problem, so that a high-resolution computation 
of the Kondo cloud with NRG is quite expensive numeri¬ 
cally. Our method remarkably deals with all length scales 
on an equal footing, so that the Kondo cloud can be 
calculated in a single step, once the many-body ground 
state is known. The coherent state expansion thus nicely 
complements the existing tools for studying the Kondo 
model. 


The rest of the Article is structured as follows. In 
Sec. [n] we define the anisotropic Kondo Hamiltonian. We 
also state the equivalent spin-boson model, and the rela¬ 
tion between the parameters of the two models. The pre¬ 
cise mapping between them is reviewed in Appendix]^ In 
Sec. |HI| we define the spatial correlation functions, that 
contain information on the entanglement cloud around 
the spin impurity, and that we will study. We first con¬ 
sider the correlations in the fermionic language of the 
Kondo model and then in the bosonic language of the 
spin-boson model. The bosonic observables are derived 
from the fermionic ones in Appendix Sec. IV intro¬ 


duces the systematic coherent state expansion. The in¬ 
novation we mentioned in the previous paragraph, that 
reduces the number of variational parameters to solve for 
numerically, is presented in Sec. |Vj In Sec.jV^we take the 
thermodynamic limit for the energy, and express it in 
terms of the remaining variational parameters. The cor¬ 
relation functions that we study are expressed in terms of 
the remaining variational parameters in Sec. |VH| Next, 
we find the optimal values for the variational parame¬ 
ters by minimizing the energy, and we present our re¬ 
sults in Sec. |VHI| We hrst confirm the accuracy of the 
method and then study in great detail the various cor¬ 
relation functions describing the Kondo cloud, as well 
as the Kondo wavefunction overlaps. Finally, in Sec. |IX| 
we summarize our main findings and, where applicable, 
compare them to results in the literature. 


II. KONDO AND SPIN-BOSON 
HAMILTONIANS 


In this section we write down the two related models 
that we study, and make explicit the deep connection 
between them. The first Hamiltonian is the anisotropic 
spin-1/2 Kondo model in one dimension, that describes 
a spin-1/2 impurity coupled to a Fermi sea. The Fermi 
sea does not necessarily represent electrons confined to 
one dimension. If the electron gas is higher dimensional, 
but the impurity is a point-like scatterer, it only inter¬ 
acts with electronic s-waves. Such s-waves in the electron 
gas can then be described by a one-dimensional Kondo 
Hamiltonian, where the spatial coordinate is the radial 
distance from the impurity. The second Hamiltonian is 
the ohmic spin-boson model, that appears in many area 
of physics, from superconducting nano-circuits to biolog¬ 
ical systems. Below an ultraviolet energy scale set by the 
Fermi energy, the two models are equivalent The 

full mapping is reviewed in Appendix [A| 

The anisotropic spin-1/2 Kondo model in one dimen¬ 
sion reads H = + H\\ + H±, where: 


Ho 


Hu 


H± 


ka 








+ 


( 2 ) 


Here Cka- annihilates an electron of wave number k and 
spin direction cr on a ring of length L, and ipa-ix) = 
J2k e.^’^^Ckal'/L. We denote these operators with tildes, 
because we will denote the useful slow modes (dehned be¬ 
low) without tildes. We assume Sk = E-fc, so that there 
are both left and right-movers at the Fermi energy. The 
superscript B indicates the bare values of the coupling 
constants and Jf. 

Integrating out high energy degrees of freedom towards 
the Fermi surface, we end up with a linear dispersion 
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relation, and an effective low energy Hamiltonian: 


Ho = 




ka 

n - '^11^ 
^11 - 

Hi_ = Ji 


|(o)y(o)-v-|(o)yyj, 
y(o)y(o)cr" +y(o)y(o)(T+ 


( 3 ) 


in units where the Fermi velocity vp = 1. We have de¬ 
fined slow modes 


Cfccr — 


Cfccr — 


(.^kp‘-\-k,a' C—fcp—/c,cr) ; 
i^kp + k,cr — C-kp-k,a) ■ 




^ikx —a\k\/2 


^kcr 7 


( 4 ) 


and a similar relation between ipa(x) and Cka- Here l/o 
is the ultraviolet cut-off energy. Operators without over¬ 
bars are associated with even-parity single particle wave 
functions, so that c\^ creates an electron in the state 
^J2/Lcos[{kF + k)x\. In D > 1 dimensions, an equiv¬ 
alent role is played by s-waves. For negative x, '0j(x) 
creates an electron in an even-parity state consisting of 
two wave packets centered around ±x and propagating 
in opposite directions towards the impurity. For positive 
X on the other hand, it creates an electron in a state con¬ 
sisting of two wave packets propagating away from the 
impurity. Because create electrons lo¬ 

calized to the same physical regions in space, we refer to 
([^ as the unfolded representation. Operators with over¬ 
bars are associated with odd-parity single particle wave 
functions, so that creates an electron in the state 
i^y2/Lsm[{kF -|- k)x]. Since the even and odd modes de¬ 
couple, and only the even modes couple to the impurity, 
we focus on electrons in even orbitals and drop the c^c 
terms in Hq. 

In the process of linearizing the spectrum (integrating 
out fast modes), the bare couplings J|^ and jf are renor¬ 
malized to new values Jy and Jl, that depend on the ul¬ 
traviolet scale I/a. Note however that SU(2) symmetry, 
if present, is preserved under renormalization. Thus, if 
jf = J|^ (isotropic limit), then J_l = Jy. 

In this work, we will make extensive use of the well- 
known fact, reviewed in Appendix that the Kondo 
Hamiltonian (|^ can be mapped onto an ohmic spin- 
boson model: 

Hsb = '^qblbg-'^^{bl + bg)s^ + (5) 

q>0 g>0 ^ 


Here hq are bosonic operators such that [6q,6g'] = 0 and 
[bq, 6^,] = Sq^qi, while Sx and are Pauli matrices. These 
pseudo-spin operators are not the physical angular mo¬ 
mentum operators of the Kondo impurity, but are related 
to them. The parameters of the spin-boson model are re¬ 
lated to those of the Kondo model by 


III. SCREENING CLOUD OBSERVABLES 


Having defined the system we investigate, we now write 
down the observables that we will study. We are inter¬ 
ested here in ground state correlation functions between 
the impurity spin and the electron spin density at a dis¬ 
tance X from the impurity, which we refer to collectively 
as the Kondo screening cloud (or cloud for short).^ We 
give equivalent expressions for the cloud correlators as 
ground state expectation values in both the fermionic 
and the bosonic pictures. In Appendix the bosonic 
expressions are derived starting from the fermionic ex¬ 
pressions. 

In the fermionic picture, the cloud correlation func¬ 
tions is defined as: 

All(x)=4(.5rP5/Ha^))j,, 

A^(x) = 4 (S^rsfix) + , (7) 

where 5*™^ = Cfjl2 is the j G {x,y,z\ component of the 
impurity spin operator, and 

Sf{x) = ( 3 ) 

(7(T' 


is the j component of the electron spin density at x. The 
subscript K indicates that the expectation value is with 
respect to the fermionic Kondo ground state. We refer to 
All as the longitudinal and A-*- as the transverse cloud. 
Both consist of a component A/ that varies slowly on 
the scale of the Fermi wavelength, and a component that 
oscillates with wave vector 2kF, and has an amplitude 
X^kp ■ The former is the result of scattering events that 
change the electron momentum by an amount that is 
small compared to kp, while the latter results from scat¬ 
tering between the Fermi points at ikp. Explicitly, one 
has: 
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X-^{x) = Xq{x) + cos{2kFx)X2k^{x), 

X'^'^{x) = ^o(a:) + cos(2kpx) — ^ sm(2kFx) X^jj^^{x). 


(9) 


For X ^ a, both the transverse and longitudinal 2kp components are proportional to cos(2kpx). In the bosonic 
language, and in the thermodynamic limit |a::| ^ L, one finds using standard bosonization identities: 


Re + (a: ^ -x), 


X^(x) = - p- 

^ ^ ’ Tx[x^ + a^) 

~ —7—S- ^ ^ ^ 

7r(a;2 + a2) 


pVH^) g-vix] 


4k = 


2x 


dm ( s. 


7r(a^+4a;2) \ ^ \a — ix 


a + ix 


gis(x) 


( 10 ) 


with the bosonic fields: 


‘P{x) = 2^ 

ij>0 



[cos(ga;) - 1] 6,, 


9{x) = 2^ 

q>0 



siTi{qx)bq. 


( 11 ) 


At a: < a, the behavior of these expressions depend 
strongly on ultraviolet physics that our model does not 
attempt to represent accurately. The regime of physical 
significance is a: a where the behavior of the cloud is 

insensitive to ultraviolet details of the model. None the 
less, we delay taking the a/x —>■ 0 limit until the end 
of the calculation, reasoning that it is of some interest to 
see how a particular cut-off scheme regularizes ultraviolet 
singularities. 

We presented expressions for a one-dimensional cloud. 
However, the generalization to a higher dimensional elec¬ 
tron gas is trivial, assuming a point impurity that there¬ 
fore only scatters s-waves.® In this case, one interprets 
the coordinate a; as a radial distance and divides the one¬ 
dimensional correlator by the area of a spherical shell of 
radius x. One also reverses the sign of the 2kp compo¬ 
nents, because the radial part of the higher dimensional 
problem is defined on the half-line, with a tt phase shift 
between left- and right movers. 

To make further progress in the computation of the 
cloud observables, we need an accurate approximation 
for the ground state of the spin-boson model. That is 
the topic of the next section. 


IV. COHERENT-STATE EXPANSION OF THE 
SPIN-BOSON GROUND STATE 


W e will use here a systematic coherent-state decompo- 
sitiorP^HSl of the many-body ground state of the Ohmic 
spin-boson model. Physically, coherent states appear as 
natural degrees of freedom, because, in the absence of 
spin-tunneling A, the two spin configurations are asso¬ 


ciated with a disp lacement of the bath modes from 0 to 
±fq = ±y/TTa/Lq. Indeed, at A = 0, there are two de¬ 
generate ground states 

I/+) O It), I/-) O It), (12) 


where 


|/±) = exp 


^'^fqibl-bq) 

. q>0 


| 0 ), 


(13) 


|0) is the bosonic vacuum, and {|t), It)} are the spin 
eigenstates of s^. The single coherent-state Ansatz (usu¬ 
ally dubbed the Silbey-Harris state) includes the effect 
of tunneling by promoting fq to a variational parameter 
and taking as trial state the linear combination |/) = 
(I/+) 0 It) ~ I/-) ® lt))/V2. This results readily in the 
simple expression = yjpaq/Lgq/{q -|- Xp) with the 
self-consistency condition Xp = X (/+|/_) = A (/| — /). 
This trial wavefunction provides an excellent approxima¬ 
tion to the true ground state provided the shifted equi¬ 
librium positions of most oscillators are not too far 
apart, i.e. if a is small compared to unity. 

For larger a, the overlap (/| —/), with fq determined 
variationally, becomes exponentially small, and this 
strongly underestimates the tunneling energy A {ux) /2. 
The coherent-state expansion addresses this issue by ex¬ 
tending the Silbey-Harris form to a more general linear 
superposition of coherent states: 


4- It)- 


(14) 
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with the set of displacement fq 
of Silbey-Harris states: 


parametrizing a family 


/(-))= exp 


Y.ftHbl-bq) 

.q>o 


|0). (15) 


Here M is the maximal number of allowed coherent states 
in the decomposition (14), and sets the level of approxi- 

and coefficients Cm are 


mation. The displacements fq 
determined by minimizing 


(m) 


£ = (V'liJsB IV') - A ('01 '0) 


(16) 


where A is a Lagrange multiplier that is used to enforce 
normalization. These parameters are generically found 
to be real in the ground state. 

The coherent-state expansion already dramatically 
improves the estimate for the tunneling energy 
A (01 10) /2 for M = 2, by allowing for cross-terms 

^/(™) I I/(’’■))^ in which the displacements and 

are anti-correlated at low q. Such cross-terms are 
therefore not exponentially small in a, and allow a siz¬ 
able energy gain compared to the Silbey-Harris approxi¬ 
mation. 

It is important to note that an arbitrary state of the 
spin-boson Hamiltonian can be written as an infinite but 
discrete sum of the form (14), at least if the bosonic bath 
contains a finite number of modes. This follows from a 
theorem, proved by Cahill,12Sl that for a single bosonic 
mode, countable sets , to = 1, 2, ...} of (real) 

coherent states exist, that form a complete basis. For 
a finite number of bosonic modes, a general state can 
therefore be approximated to any required accuracy, by 
making M sufficiently large. It is in this sense that the 
coherent state expansion is numerically exact, provided 
in practice that good convergence to the true many-body 
ground state occurs. Previous investigation^^, as well 
as the extensive comparisons made in Sec. |VHI[ demon¬ 
strate indeed that the expansion (|I4[) rapidly approaches 


the exact ground state of the spin-boson model as M is 
increased to moderate values, also for an infinite bath. 


V. REDUCING THE NUMBER OF 
VARIATIONAL PARAMETERS 

In previous implementations of the coherent expan¬ 
sion ( I4p5 ), the bosonic bath was restricted to a large 
but finite number N of modes, and each of the M x N 
displacements was treated as an independent vari¬ 
ational parameter. In this work, we want to calculate 
the Kondo screening cloud, which requires a high spatial 
resolution from short to possibly exponentially large dis¬ 
tances, so that a huge number of bath modes needs to 
be included. It is therefore desirable to have an imple¬ 
mentation of the method in which the number of vari¬ 
ational parameters does not depend on the number of 
included modes. Such an implementation would also al¬ 
low for computations for a continuous spectral density 
(thermodynamic limit), which numerical techniques like 
NRG are not able to perform. In this section, we derive 
such an implementation, which is the the main technical 
innovation of our work. 


We show that the algebraic structure of the variational 
equations constrains the functional dependence of the 
displacements with respect to momentum q to such 
an extent that there are in fact only -|- M — 1 varia¬ 
tional parameters, independently of the number of modes 
N of the bosonic bath. To stress that the argument does 
not require a linear bath spectrum, we replace the term 
^g>o qb^qbq in the spin-boson Hamiltonian, with the more 
general kinetic energy expression J2q>o 
mainder of this section. In subsequent sections, we will 
specialize again to the linear bath spectrum. 

Exploiting the coherent state structure of the Ansatz, 
we can express the total energy £ in terms of /g"^^ and 
Cm, 3'S. 


M r 

£= 'E (/(-)| [a;g/M4" 

m,n—l \ g>0 


9 q 

2 


(4™^+4”^) - f > (17) 


where 


with the use of 0 and ([T^, this leads to 


j(m) 


±/(”)) = exp 




q>0 


■ (18) 


M 


^ {U^uUJq + 4”^ = Wm9q, (19) 


n—1 


Considering the minimization condition d£ jdf^E — 9 


where the entries of the M x M matrices U and U, and 
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the M-dimensional column vector W are: 

Ulm = 2Cm 

^Im — Plm P Qlm ^Im ^ ^ {_Pln Qln') ; 

n 

X E - f (/t> -/r) - ^ 

k>0 


Qlm = 


( 20 ) 


in considering a bath with a very large number of modes 
has been removed. The next question is whether the en¬ 
ergy, which involves sums over all bath modes, can be 
calculated efficiently for a large number of modes. Only 
single sums over the momentum have to be computed in 
the energy functional, and in the worst case, the total 
numerical cost is linear in the number of modes. This is 
a major improvement with respect to a brute force di- 
agonalization of the model, where scaling of the Hilbert 
space dimension is exponential with the number of de¬ 
grees of freedom. In order to reach the continuum limit, 
the discrete momentum sums are replaced by integrals, 
which fortunately can be performed analytically rather 
than numerically, as we now demonstrate. 


It is crucial to note that U, V and W in (19) do 
not depend on the value of the momentum q labelling 
the displacement that we consider within the varia¬ 
tional equation. These matrices and vectors however de¬ 
pend non-linearly on the complete set of displacements, 
through sum mations ove r a dummy momentum index as 
in Eq. (18). From (19), it readily follows that is 
given by: 


M 

+ ( 21 ) 

n—1 

Since Uu}q-\-V is linear in Wq, 

]\f{m) ( \ 

J2HUu,, + Vr‘]„,„W,. = ^^, ( 22 ) 

where is a polynomial of order M — 1 in z and 

D{z) = det(t/z -I- H) is a polynomial of order M in z. 
Note that the denominator D{z) is common to all sets 
of displacements for m = 1,... ,M. Using the ex¬ 
pression ( [^ for W in terms of U, we see that for large 

z, —>■ l/2z independent of m. Thus we arrive at 

the form: 


fir 

Jq 

hjyj 




M-1 


■s^M - 

(?) = 




(23) 


where /imM-i = 1, and in general /imn is real. The poles 
uin that are not real come in complex conjugate pairs. 
The optimal state I'i/') can be obtained by numerically 
minimizing the energy with respect to the M x {M — 1) 
unknown coefficients the M poles a;„, and the M—1 
weights Cm (accounting for wavefunction normalization), 
so that there are A'P + M — 1 parameters to be found in 
total, independently of the number of bosonic modes in 
the problem. 

Since the dimension of the search space does not de¬ 
pend on the number of bath modes, the main obstacle 


VI. THERMODYNAMIC LIMIT AND THE 
ANALYTICAL EVALUATION OF MOMENTUM 
INTEGRALS 


In this section we specialize again to the case of a linear 
bath spectrum Wq = q, and analytically perform the mo¬ 
mentum integrals involved in the calculating the energy 
for given set of variational parameters ^m and Cm- 
We do so for the smooth ultraviolet cut-off gq oc 
that the bosonization of the Kondo model naturally intro¬ 
duces. Although we do not work out the details here, we 
note that the momentum integrals can also be performed 
analytically in the case of a sharp cut-off gq cx 9(l/a — q). 
In the L —>■ oo (thermodynamic) limit, momentum sums 
are replaced by integrals according to 


27r 

T 


E- 

q>0 



(24) 


With fq ' of the form (23), and gq of the form (|^ dic¬ 
tated by the mapping from the Kondo model, the energy 


(•^l ffsB \i’) 

(V'l V') 


(25) 


can be written as a functional of the variational param¬ 
eters as: 

E"„, - Ae-J'i;) 

^ =- zzjiZ. - -• ™ 

which involves the Laplace transforms 

pOO 

(27) 

JQ 

of the rational functions 

J^lik) = k[hm{k)±K{k)]\ 

J^mJk) = [khm{k)-l][kK{k)-l]. (28) 

Note that we have omitted here a constant term 

= ( 29 ) 















Im k 


for X > 0, and 



Re k 


FIG. 1: Chosen deformation of the contour for the inte¬ 


gral (311. 


that does not depend on the variational parameters. The 
functions have M second order poles of the type 

[k — . For large k, the behavior is 


^ k 


-2+A 


(30) 


The Laplace transforms can be performed analytically, 
using the following method. In general, we consider in¬ 
tegrals of the form 


I = 


pOO 

/ dke-^'^R{k), (31) 

Jo 


where the function R{k) is a ratio of polynomials, and all 
poles are of second order, i.e. 


Rik) = 


p{k) 




(32) 


This excludes the possibility of poles on the positive real 
line, which yields E = -|-oo. Furthermore, the numerator 
p{k) is a polynomial of at most order 2M — 1. 

It would be useful if we could replace the integration 
contour in (311 by a closed contour such as the one in Fig¬ 
ure We can do so, provided we multiply the integrand 
by a function J-{k) with a branch cut along the posi¬ 
tive real line such that E{k + t0+) — E{k — i0+) oc 9{k). 
Furthermore, J- must be such that the contribution to 
the integral that comes from closing the contour at large 
k is negligible. A function that exactly meets these re¬ 
quirements is r(0, —ak), where r(0,z) is the incomplete 
Gamma function of order zero.^ It is defined as 


r(o,z)= r 

J z 


dt—, 


(33) 


where the integration path does not intersect the nega¬ 
tive real line and excludes the origin. Its main relevant 
properties are 

r(0, -(x -b i0+)) - r(0, -(x - i0+)) = 27rf, (34) 


r(o,-z)~--, 

Z 


(35) 


for large \z\. Furthermore r(0, —z) is analytical except 
in the neighborhood of the positive real line. Thus we 
find 

r°° f flk 

/ dke-^'^R{k)= / -i?(fc)e-“'=F(0,-afc) 

Jo Jc 27ri 

= ^Res{i?(A:)e"“''r(0,-a/c),a;„} . (36) 

n 

In the first line, C is the contour depicted in Figure 
Returning to the integrals , one then has 


M 


= ^Res{jA„(fc)e “'=r(0,-a/c),a;z} . 


/ = 1 


Setting 


Fiz)=e-^T{0,-z), 


and noting that 


= - 


r(j) + - 


evaluation of the residues yield 

iLu = ^mn/^(wza) - ( aF{uJia) 

i 

where 


OJl 


(37) 

(38) 

(39) 

(40) 


Kmni = lim(fc-a;i) J„„(fc), 

k—^uji 

Kinnl = ^ [{k - Jtnik)] ■ 


(41) 


To relate the K parameters to the coherent state dis¬ 
placements explicitly, we define 


hmn = \lTSl {k - UJn)hm{k) 
k^UJn 

^M-1 


E m—1 I 

1=0 dml^n 


(42) 


and 


hran = Hm — [{k - UJn)hm[ki)] 
k—GiK 


E m — 1 J 
I—I 


1-1 


- K 


^ E 

l — l^n 


UJn - 


(43) 


This gives 


K: 


±1 

mnl 


= OJlQlml F hniY 


k^mnl — ^^ihjnlhnh 


^0 

‘■rr^ 

AO 


k^mnl — (kmlFhnl) + 2uJl{hml F hnl){hml F hnl), 


k^mnl UJih^i l)hyj/ 

+ ^ihmlihnl F UJlhnl ~ !)• 


(44) 
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Substitution of (441 into (|40[), and (401 into (26) yields 


an explicit expression for the energy in terms of the vari¬ 
ational parameters. 

As a consistency check, we can work out the single co¬ 
herent state theory (M = 1) to make sure that we recover 
the results of Silbey and Harris. This approximation is 
accurate for a <C 1, which, in the Kondo language, corre¬ 
sponds to the strongly anisotropic Jy ^ Jl regime. For 
M = 1, there is only one variational parameter, which 
is conventionally defined as A/j = —wi. This tra nslates 
to hii = 1 and An = 0 [see Eqs. (42) and (43)]. The 
energy of the Silbey-Harris state then reads 


E = -- [aA?j9A«F(-Ai^a) -y A (/| -/)] , (45) 


with 

(/I -/) = exp {-a [F(-Aija) -f Ai^aAR-F(-Ai^a)]} , 

( 46 ) 

and we have used the identity ( [^ to relate F and its 
derivative. Minimizing E with respect to An, we recover 
the known self-consistency condition 

Afi = A(/|-/). (47) 

Assuming than An ^ 1/a, and using the fact that for 
small z, 


F{-z) = -\n{z)--fE + 0{z), (48) 

where je = 0.577 is the Euler-Mascheroni constant, one 
can solve the self-consistency condition for An to obtain 

An = A{e^+^^aA)^ . (49) 

Up to a pre-factor. An corresponds to the Kondo temper¬ 
ature in the regime of small a, where the single coherent 
state Ansatz is accurate. 


In general, it is far more efficient to evaluate the energy 
analytically, as was done in this section, than to do the 
integrals in (26) numerically. There are however small re¬ 
gions of the search space where the analytical evaluation 
of the energy is not stable. These are regions in which 
two or more of the poles tOm lie close to each other. In 
these regions, there are large cancellations between indi¬ 
vidual terms in the sum over residues in (37), and this 
leads to large numerical errors when individual residues 
are calculated separately before they are summed. We 
circumvent the problem as follows. When the minimiza¬ 
tion algorithm searches a dangerous region of the search 
space, it does not try to evaluate the residues at the of¬ 
fending poles individually. It rather takes the slow but 
safe option of numerically integrating around a loop that 
circles all closely spaced poles at a safe distance. Fortu¬ 
nately, one does not have to fall back on this contingency 
plan too often, as the problematic regions of the search 
space are small and do not seem to be particularly fa¬ 
vored in the actual optimal solution. 

VII. EXPRESSING THE CLOUD WITH 
COHERENT STATES 


In the previous section, we obtained an analytical ex¬ 
pression for the energy in terms of M‘^ + M — 1 variational 
parameters. This result allows for a significant speed-up 
of the numerical minimization of the energy. In this sec¬ 
tion, we apply the same analytical technique to evaluate 
the momentum integrals involved in the calculation of 
the Kondo screening cloud, in terms of bosonic displace¬ 
ments. Evaluating the four different cloud correlators 
(10) for the M-coherent state wavefunction can be done 
using straight-forward coherent-state algebra: 


M 


Xoix) = - 


cos(yaH(a;)„„)cmCn (/ml-/n) 




7r(x^ + a^) 

(x) = — , (/ml -/„) 

7r(x^ + a^) 

' ' m n — 1 


m,n—l 

M 


m,n=l 

a 


M 


= -(1 - V^) / 2i 2'' + 9 C{x)mnCmCn{f 7 n\ fn) , 

7r( 2/ (X ) ZilT 


2x 


X" fri = 

7r(4x^ -I- a^) 


Im 


a + ix 


a — IX 


M 


Y. e*^^(")-CmC„(/m| /n) 


together with the normalization condition 

M 

^ ^ ^mC-n (/m| fn) — !■ 
m,n=l 


(50) 


The matrices A, B, C, and D are defined as: 

pOO 

(51) A(x)mn = J dq cos{qx)e~°‘‘^[hm(,q) + hn{q)], 

pOO 

B{x)mn = / dq s\n{qx)e~'^'^[hm{q) - ^n(g')], 

poo 

C{x)mn = dq cos{qx)e~'^'^[qhm{q) + qhniq) - 2], 
^0 
pOO 

D{x)mn = dq sin{qx)e~'"^[hm{q) + hn{q)]. (52) 

Jo 
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Note that B, C and D remain finite if the limit a —>■ 0 
followed by a: —>■ 0 is taken, while A diverges logarith¬ 
mically. The method we used to evaluate {fm\ ±/n) = 
exp(—a/;Jj^/4) analytically can be extended to evaluate 
the above integrals as well. The detail of the calculation 
can be found in Appendix The resulting expressions 
are: 


M 

Ai^x') jYin — Re ^ ^nl)B i^UJi (u ; 

l=l 
M 

B{x)mn = “ hni)F{uJi{a - ix)), 

1^1 

M 

C{x)mn = Re'^u;i{hmi + hni)F{uJi{a - ix)), 

1^1 

M 

D{x)mn = + hni)F{uji{a - ix)), (53) 

1^1 


where F{z) is defined in (38). 

For future reference, we now consider the large (com¬ 
pared to the Kondo length) x asymptotic behavior of the 
expressions (50) for the screening cloud. At \z\ » 1, 


- 4 - 


(54) 


The residue theorem can be used to derive the identities 

M M 










(55) 


With the aid of the above equations, the cloud correlators 
(50) are found to decay like Xjx^ at large x. Explicitly, 
one finds for x ^ : 


= X^{x) = 


X",^{x)=X"{x) . 


X^ 

a _ C|| 


X^ 


(56) 


where the lengths ^_l and are given by 

M 

f, — Z_ W f,V^I2Zl(f^rr.t+hr,l)Flia) Ir \_f\ 

S-L — ^ ^ <-m^n \Jm\ Jnj 5 

m,n=l 
!— M 

^\\ = F F, + hniO)]CmCn {fm\ fn) ■ (57) 


The above expression establishes a very direct connec¬ 
tion between the q = 0 behavior of the coherent state 
displacements and the large x behavior of the longitudi¬ 
nal cloud. The small contribution Oiajx^) to Ag is a 
feature of the smooth ultraviolet cut-off ~ 6““"^, and will 
be ignored in what follows. Note that the Okp and 2kF 
components of the cloud become equal at large x. Fur¬ 
thermore, there are two emergent length scales ^_l and 
^11 in the problem, which are in general not equal, unless 
spin-isotropy is restored. 

In the single coherent-state approximation, which is 
accurate for small a, the cloud correlators can be explic¬ 
itly computed: 


X^{x) = a4^(x) = - 
= -(1 - Va) 


aAn 


2v^[F(-AHa)-Re F{-An(a-ix))] 


fA{x‘^ + a^) 
a ^faAR 


f(xA + a^) TT 


Re F'(—Ai?(a — ix)), 


= 


2x 


7r(a2 -I- 4a;2) 


Im 


O' F IX ^2iv^Im F{-Aft{a-ix)) 


(58) 


where, as noted in Sec. Ar = —uji is the Kondo 
energy scale. The two correlation lengths and [cf. 
(57)] are in this anisotropic limit: 


_ oAr 2^F(-ARa) 
^^~2fA 


C|| = 


ttA 


R 


(59) 


These simple analytical expressions for the Kondo cloud, 
valid in the limit of strong spin-anisotropy, have not ap¬ 
peared in the literature before. We will analyze them 
further in the next section, along with numerical results 
obtained for larger a values from the systematic coherent 
state expansion. 


VIII. RESULTS 


In Sec. |VI|and Sec. m we have collected the tools to 
calculate the average energy and the components of the 
screening cloud, for an M coherent-state wavefunction, 
in terms of only 0{M^) variational parameters. In order 
to find the variational parameters, the energy must be 
minimized, and for this p urpos e we use a standard simu¬ 
lated annealing algorithm .ESEU Clearly, the quality of the 
approximation is limited by the maximum number of co¬ 
herent states that can be handled with the available com¬ 
putational resources. Using a single personal computer 
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FIG. 2: (color online) Fractional improvement {Em — 
Em-i)/Em -1 in the minimum energy Em as the number 
of coherent states M is increased, at three different points 
in parameter space. Circles correspond to a strongly spin- 
anisotropic point (a = 0.3, A = 0.05/a), triangles to a mod¬ 
erately spin-anisotropic point (q = 0.6, A = 0.05/a), and 
diamonds to a perfectly spin-isotropic point (q = 0.85, A = 
0.156/a), associated to equal Kondo exchange couplings J± = 
J|| = 0.49. 


and simulated annealing minimization, we have found it 
possible to go up to M = 7 coherent states, which is 
enough for our purposes here, although an improved al¬ 
gorithm combining global and local optimizatiorP^ can 
reach values as large as M = 24. In fact, the required 
number of numerical operations is not a limitation per 
se here. Rather, the main difficulty is that the energy 
landscape in the space of variational parameters is very 
shallow and contains several low lying minima. 

This section is divided into three extended sub-parts. 
First, we benchmark the coherent state method, both 
against previous results, and also by studying the con¬ 
vergence properties of the coherent state expansion, es¬ 
tablishing its domain of validity for the computation of 
the screening cloud. Then, we consider the Kondo over¬ 
laps proposed in Ref. |5T1 which we compute both ana¬ 
lytically in the spin-anisotropic regime, and numerically 
for the nearly spin-isotropic case. Our results quantita¬ 
tively agree with recent field theoretic calculations.!^ We 
finally perform an extensive study of the Kondo cloud 
correlators over a wide range of parameter values by a 
combination of analytical arguments and extensive co¬ 
herent state simulations. This allows us to uncover the 
precise universal features of the screening cloud of the 
anisotropic Kondo model. 


A. Convergence properties of the coherent state 
expansion 


FIG. 3: (color online) Energy uncertainty AE as a function 
of the number M of coherent states, for the same three points 
in parameter space as the data in Figure]^ 


M of coherent states increases, using various observable 
quantities. In Refs. and [531 convergence was estab¬ 
lished only for a spin-boson model containing a large but 
finite number of modes, and we demonstrate here that 
fast convergence also occurs for a continuous bath spec¬ 
trum. 

We start by plotting in Figure the fractional im¬ 
provement {Em — Em-i)/Em -1 in the minimum energy 
Em obtained by adding an extra coherent state in the 
Ansatz, as a function of M. Note that, as in Sec. VI we 


do not include the constant Eq = —aj^a in the defini¬ 
tion of Em- This convention leads to a denominator in 
{Em —Em-i)/Em -1 that is closer to zero, and therefore 
to a more stringent measure of convergence. Figure 
shows results for three points in parameter space. One 
of the points, {a = 0.3, A = 0.05/a), corresponds to a 
strongly anisotropic situation where convergence is very 
rapid. Another curve {a = 0.6, A = 0.05/a) corresponds 
to a less anisotropic situation, where the convergence is 
slower. The remaining point {a = 0.85, A = 0.156/a) 
corresponds to the isotropic coupling J_l = Jy = 0.49, 
or a Kondo temperature Tk = exp(—7r/J)/a ~ 1.6 x 
10“^/a. For (a = 0.3, A = 0.05/a), the minimum en¬ 
ergy changes by an amount comparable to the accuracy 
goal of the minimization module, by the time that M 
reaches 6. For the other two points in parameter space, 
the marginal change in the minimum energy is a fraction 
of a percent at M = 7. All the results that we present 
below are at least as converged as these last two cases. 

As a complementary test, we can also verify that the 
multiple coherent state \'i/i) converges to an eigenstate of 
E[ as the number of coherent states increases. To do so, 
we calculate the energy uncertainty 


AE = 


I IV’) - 


\H\ 


(60) 


In this subsection, we present strong evidence for the 
rapid convergence of the wavefunction (14) as the number 


for the optimal coherent states Ansatz l"^), as a function 
of M. For this purpose, an expression for the overlap 
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FIG. 4: (color online) The Okp transverse component X(^ of 
the screening cloud, versus distance x from the impurity, for 
a = 0.1, A = 0.05/a. The black dashed curve is the (fully 
converged) result of a 6 coherent state calculation. The solid 
blue curve corresponds to the single coherent state result (581. 
In the spin-anisotropic limit a <C 1, the single coherent state 
approximation is thus very accurate. 


FIG. 5: (color online) Convergence of the cloud for the Okp 
transverse component versus distance x, for a. = 0.85, 
A = 0.156/a. This parameter choice corresponds to an 
isotropic Kondo coupling Jx = J\\ = 0.49. The black dashed 
line represents the result of a M = 7 coherent state calcula¬ 
tion, while solid curves represent results for M = 2, 3 and 5 
coherent states. 


^/(m)| ^2 |j(n)^ ig needed, 
integrals as in Sect. VI are involved. 


It turns out that the same 
In terms of 


these integrals, the overlap is expressed analytically as: 





+ (7+1 -/-M 

^ g '-'a \^mn ^mn) ' 


(61) 


The rest of the calculation can then proceed with the 


technology developed in Sect. VI for evaluating the in¬ 
tegrals Note that derivatives dal^m refer to par¬ 

tial derivatives with respect to the explicit a dependence 
of the integrals, and not to total derivatives involving 
the implicit a-dependence of the variational parameters. 
The energy uncertainty XE is showed in Figure as a 
function of coherent state number M, for the same three 
points in parameter space for which we have investigated 
the energy convergence above. In all three cases the en¬ 
ergy uncertainty decreases monotonically, and roughly 
as a power law oc The fact that the uncertainty 

clearly tends to zero as M increases, shows that the trial 
wavefunction © converges to the true ground state as 
the number of coherent states is increased. 

Next, we investigate the convergence of the screening 
cloud itself. In Figure we show an instance of very 
rapid convergence in the case of strong spin-anisotropy, 
and in Figure]^ an example where there is a noticeable 
change between M = 1 and M = 7 coherent states (note 
also that this data set is one of the least converged ones 
included in this work). In both figures we plot the Okp 
transverse correlator Xq{x), as we found this component 
to show the most dramatic change as M is increased. 
Figure 1^ corresponds to a = O.I and A = 0.05 (strong 


spin-anisotropy), and the comparison to a M = 6 co¬ 
herent state calculation shows that the single coherent 
state (Silbey-Harris Ansatz) is nearly exact in this limit. 
The results in Figure were obtained at a = 0.85 and 
A = O.I56/a, which corresponds to an isotropic Kondo 
coupling Jj_ = J|| = 0.49, for M = 2, 3, 5 and 7 coher¬ 
ent states. Good convergence is clearly ensured by the 
computation with M = 7 coherent states. 


B. Comparison to exact results based on 
integrability: ground state energy and Kondo 
overlaps 

In this subsection, we give further evidence that 
ground state properties of the Kondo model can be calcu¬ 
lated accurately using the coherent state expansion. In 
particular, we focus here on several physical quantities 
that can be computed exactly via the Bethe Ansatz or 
related integrability techniques. 

In Figure we compare the multiple coherent state 
estimate for the ground state energy to results obtained 
via the exact Bethe AnsatJ^, at A = 0.05/a and vari¬ 
ous a values. The high energy cutoff 1/a corresponds to 
the parameter Wc in Ref. 1421 and we used from this ref¬ 
erence expression (C.9) for the ground state energy and 
equation (8) for the Kondo temperature Tp, with the re¬ 
lationship between D and oJc given in (C.8). Although 
our calculation is variational, the energy is typically con¬ 
verged to about 0.1% (see Figure [^. One has to bear 
in mind however that, in the ultraviolet, the model for 
which the Bethe Ansatz yields the exact solution differs 
from the model we consider here. The Bethe Ansatz re¬ 
sult is only valid if all relevant energy scales in the prob¬ 
lem are much smaller than the ultraviolet cut-off scale. 
For finite Aa, the Bethe Ansatz expressions even present 
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a 


FIG. 6: (color online) Ground state energy versus a as ob¬ 
tained via the multiple coherent state approximation (dots) 
and the exact Bethe Ansatz solution (solid line), for A = 
0.05/a. 



FIG. 7: (color online) Kondo overlap {‘tp{J±)\ between 

ground states at different and varying values of Jx, but with 
the same fixed value of Jy, versus the ratio Tk{J_l)/Tk{J±)- 
Here, the strongly spin-anisotropic regime is considered with 
a = 0.1, i.e. J|| = 4.30. The solid line shows the single 
coherent state formula (651, whereas the dashed line is the 
exact analytical result ( |62[|. The dotted line is the small a, 
small \z\ approximation (|66||. 


spurious divergences around a = (2n -|- l)/(2n -|- 2), for 
u = 1, 2, .... As a result, exact agreement can only be 
expected at \Ea\ 1, and this explains why our vari¬ 
ational result is slightly off the Bethe Ansatz result for 
small a, although the variational approach presents bet¬ 
ter convergence in this regime. In general, the numerical 
data agrees so closely with the Bethe Ansatz result that it 
is impossible to distinguish between errors resulting from 
truncation of the coherent state expansion at M = 7 and 
errors due to ultraviolet differences between the models. 

Another interesting quantity for which an exact analyt¬ 
ical result has recently been obtaineci^, and which gives 
some indirect information on the Kondo cloud, is the 



FIG. 8: (color online) Kondo overlap ('i/'(Jx)| at in¬ 

termediate a values. Open symbols in various colors show the 
M — 6 coherent state results for Jy = 0.905 (i.e. a = 0.3) 
for several values of J± and G [0.016,0.346] (i.e. A £ 
[O.OO 5 / 0 , 0.110/a]). Closed symbols show the M = 7 coherent 
state results for Jy = 0.451 (i.e. a = 0.6) for several values of 
J± and Ji € [0.016,0.346] (i.e. A € [0.236/a, 0.785/a]). Solid 
curves represent the exact analytical result (62l. The Kondo 
temperature was taken from the approximate identification 

rx = i/Ci|. 



FIG. 9: (color online) Kondo overlap {ijj{J±)\ at a = 

0.8 (J|| = 0.663). Open symbols in various colors show the 
M = 7 coherent state results for several values of Jx and 
Ji G [0.377,0.880] (i.e. A G [0.12/a, 0.28/a]). Solid curves 
represent the exact analytical result (62l. In the main panel. 


the improved renormalization group estimate ( |64[ ) was used 
for Tk- In the inset, the less accurate value Tk = l/Cy was 
used. 


overlap of two Kondo ground states with 

different transverse exchange couplings, as a function 
of the Kondo temperature ratio T/y(Jx)/TR-(Ji). Here 
\i^{J±)) and \ijj{J±)) denote the full many-body ground 
states obtained at different J± but for the same Jy, while 
Tk{J±) and are the associated Kondo tempera- 
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tures. With the definition z = ln[TifJ_ l)], the 
exact result reads: 

/,/TM ifT'w 1 sinh[(l - q;)z/2] 

»(^x)l iiW) = 1_„ ,i„i,(,/2) »-(--)■ 


ga{z) = exp 


r°° dt sin^zt/T:) sinh(^) 
Iq t sinh(2t) cosh(t) ( t 

V 1 —Q! 


(62) 


When comparing our numerical results to this analyti¬ 
cal formula, a subtle issue arises. The Kondo tempera¬ 
ture is certainly related to the inverse of the size of the 
screening cloud, but the exact relation may well involve 
an 0(1) factor that is J± dependent. We have therefore 
tried various definitions of the Kondo temperature. For 
small to moderate a, we find that Tk oc 1/C||i as defined 
in Eq. (56), with a J_L-independent proportionality con¬ 
stant, works well. Indeed, in the Silbey-Harris regime, 
this correspondence is exact (see below). For larger a 
however, this definition seems to incur a systematic er¬ 
ror. For J|| and J_l relatively small, as is exemplified by 
the data we collected at a = 0.8, we have found it bet¬ 
ter to estimate the Kondo scale from direct integration 
of the standard weak-coupling in J± and Jy poor man’s 
scaling equations: 


dJ± 

hr 


= ^ Jii J 

TT 




^■^11 

di 


= -J 


_L5 


(63) 


which lead to an expression of the Kondo scale for the 
spin-anisotropic Kondo model that is valid for Jy ^ 1 
and J± <C 1: 


Tk = - exp 
a 


Jl - 4 


- arctan 


IJl-jf 


(64) 

Let us first consider the regime of small a. Here we 
have seen that single coherent state results are already 
well-converged. In the single coherent state approxima¬ 
tion, the overlap is given by: 


Ref. HU in which ln(^/i(Jj_)| is calculated to sec¬ 

ond order in the impurity interaction, and then exponen¬ 
tiated. In Figure we compare the single coherent state 
approximation to the exact result (62). we see that, un¬ 
like the semiclassical formula 


the single coherent 
state approximation remains valid up to large ratios of 
the Kondo temperatures, because the ground state is well 
captured for any value of the Kondo temperature (pro¬ 
vided a is small enough). 

We now turn to the regime of larger dissipation, and in 
Figure]^ we show multiple coherent state resnlts for the 
Kondo overlap at a = 0.3 and a = 0.6, again using the 
(now approximate) identification Tk oc I/C||) together 
with the exact result (62). Data points of the same color 
were obtained by keeping J_l fixed and varying J(|_ < Jj_. 
If two data points have different colors, they correspond 
to distinct Jl and J'j^. It is therefore already nontrivial 
that the differently colored points fall on the same curve, 
confirming the universality predicted by the exact analyt¬ 
ical expression (62). Generally, good agreement with the 
analytical result (621 is seen, with only a small system¬ 


atic error pushing the numerical curves slightly below the 
analytical ones, which increases as a increases. A small 
part of the error is likely due to an error of a few per¬ 
cent in the numerical estimation of ^y, and to Tk being 
only an order of magnitude or so less than the ultraviolet 
scale 1/a at large J±. These errors would be reduced if 
we could use more coherent states and smaller J±. How¬ 
ever, as pointed out earlier, the main discrepancy comes 
from the chosen definition of the Kondo temperature in 
Ref.EU which is not exactly equivalent to 1/^y at larger 
a. Indeed, Figure shows the multiple coherent state 
results for the overlap at a = 0.8, using both the naive 
identification Tk oc 1/^y (inset) and the renormalization 
improved estimate ( [m| ). The latter choice leads clearly 
to a substantial reduction of the error. Our conclusion 
is that the actual Kondo overlap ('0 (J_l)| 'ip{J'j_)) is cal¬ 
culated quite accurately in the coherent state expansion, 
and that the (modest) observed errors are associated with 
the extraction procedure of the Kondo temperature from 
the size of the cloud. 


(^/>(Jj_)| V’('/i)) =exp| l-|coth(z/2) 


(65) 


where z = In [Aij(J(L)/Afl( J_l)]. For small a, Ref. HT] 
quotes the “semi-classical” result. 


C. Comparison to the exact longitudinal 
forward-component of the cloud at the Toulouse 
point 


(V’(J±)IV’(Ji)) = l + f [l-|coth(z/2) 


( 66 ) 


which can be obtained by expanding the exact result (62) 


to first order in a. Since zcoth(z/2) grows linearly in z 
for large z, the semi-classical result can at best only be 
valid for z sufficiently smaller than l/a. Referring back 
to the single coherent state result (59) for the correla¬ 


tion lengths, we see that if we make the identification 
Tk oc 1/C||i the single coherent state approximation ([6 m 
is nothing but a resummation of the small a resulx in 


In this subsection provide strong evidence that the 
coherent state expansion (14) captures the full spatial 
structure of the Kondo cloud, using an exact analyti¬ 
cal result at the so-called Toulouse point. Indeed, at 
Jy = 27r(I — I/-\/2) or equivalently, a = 1/2, the Kondo 
model is equivalent to a fermionic non-interacting reso¬ 
nant level model [43] . An exact result for the longitudinal 
Okp component of the screening cloud can be obtained 
following the route set out in Refs. HUandHU We review 
this calculation in Appendix]^ The result is given by the 


























15 



FIG. 10: (color online) The longitudinal Qkp component of 
the screening cloud at the Toulouse point a = 1/2. Five 
different values of A ranging from 0.02/a to 0.1/o were used. 
The symbols correspond to results obtained with the coherent 
state expansion (truncated at M = 7 terms) for a choice of 
discrete positions, while the solid lines represent the exact 
result ( |67[ ). Also indicated for each curve is the associated 
Kondo length ^||, ranging from 18.2o at A = 0.1/a to 456o at 
A = 0.02/a. 


simple formula 

Xll(x) = 


y/2 Jj J? 1x1 


( _ 

TT^ 47ra 


47ra 


(67) 


with F as defined in Eq. (38). From the point of view 


of the coherent state expansion, there is nothing special 
about the Toulouse point. In the coherent state expan¬ 
sion, Xq is expressed as a linear combination of F func¬ 
tions with different position dependent arguments, cf. 
Eq. (53), and an infinite number of te rms are required 
to approximate F^ exactly. In Figure 10 we compare 


the exact expression to results obtained with the coher¬ 
ent state expansion, truncated at M = 7 terms. For 
clarity, we show here results of the coherent state ex¬ 
pansion only for a discrete set of x values, because the 
complete curves would have completely covered the exact 
result. We clearly find near perfect agreement between 
the coherent state expansion and the exact result. Since 
the coherent state expansion does not exploit any special 
features of Toulouse point, we expect similar accuracy at 
a similar cost (M = 7) to be achievable at other values 
of a. 


D. Detailed analysis of the screening cloud 

Having thoroughly established the accuracy of the co¬ 
herent state approximation for the Kondo ground state, 
we now proceed to investigate the physical features of the 
Kondo screening cloud, both for isotropic and anisotropic 
regimes. 


1. Summary of known results 

We first briefly recall the available analytical results 
regarding the isotropic screening cloud. In the isotropic 
case (J_L = T||), the transverse components of the screen¬ 
ing cloud equal twice the longitudinal components, i.e. 

xjl=X^/2 = Xo, xl^=X^,j2 = X2k,. (68) 

For X ^ a, both Xq{x) and X 2 kp(x) are expected to be 
universal scaling functions 

^Xk{x) = Xk{x), x = x/^, k = 0,2kF. (69) 

Here Xk is independent of J, and all parameter depen¬ 
dence is contained in the Kondo length The Kondo 
length is expected to be inversely proportional to the 
Kondo temperature, but the exact relation may contain 
a J-dependent proportionality factor of order unity, as 
we discussed previously. The following asymptotic re¬ 
sults for Xq and X 2 kF have been derived analyticallj® 

Xo{x) oc - , , X 2 kp{x) oc ' / for i < 1 

xymx)^ xlnx 

Xo(x) ~ X 2 /cj,(i) ~fori>l. (70) 

The regime of small x is perturbatively accessible with 
a calculation using renormalization group techniques, 
while the large x regime is treated using Fermi liq¬ 
uid theory. Note that at small x, the 2A:F-oscillatory 
component of the cloud dominates the Okp component 
slightly. This implies that the total correlation function 
X{x) = Xo{x) + cos{2kFx)X2kF oscillates between posi¬ 
tive and negative values, with wavelength -K/kp at small 
x. In other words, close to the impurity, the spin cor¬ 
relations between the electron gas and the impurity al¬ 
ternate between being ferromagnetic and being antiferro¬ 
magnetic, on the scale Tr/kp- The crossover from slower 
than 1/x decay inside the cloud to 1/i^ decay at large 
X is expected to occur at x ^ i.e. x ^ 1. Our goal 
in the next paragraph is to confirm these results for the 
isotropic cloud, and to compute the full crossover curve 
from the coherent state expansion. In a second step, we 
will examine how the cloud correlations change when the 
Kondo couplings are not isotropic (J_l ^ Jy). 


2. Kondo cloud in the isotropic case 

As a first step, we present results in Figure [TT] for the 
screening cloud calculated for several parameters on the 
isotropic line J± = Jy, or equivalently a = (1 — aA/2)^. 
The correlation functions are plotted in units of 1/a on 
the vertical axis and in units of a on the horizontal axis. 
When considering the small x regime numerically, one 
must remember that x has to remain sufficiently larger 
than the short distance ultraviolet cut-off a, since the 
limits a —>■ 0 and x —> 0 do not commute. This restriction 
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FIG. 11: The four components of the isotropic screening 
cloud: the top panel represents the longitudinal (black) and 
transverse (gray) Qkp components, while the bottom panel 
represents the longitudinal (black) and transverse (gray) 2k f 
components. The calculation was performed for the choice of 
Kondo couplings J — 0.42, 0.49, 0.66, and 0.84, which corre¬ 
spond to a = 0.87, 0.85, 0.8, and 0.75. The various curves can 
be identified from the fact that the crossover to faster 1/x^ 
decay occurs on an increasing length scale ^ as q is increased 
(or equivalently J is decreased). For a = 0.75 convergence 
was achieved with M = 5 coherent states, while for the other 
values of a, M = 7 was required. 


is implied whenever we consider the asymptotic a; —> 0 
limit. The transverse and longitudinal components are 
plotted in the same panel, and we expect x\ = j2 
for both the k = Okp and k = 2kp components due to 
strict spin-isotropy. The numerical results of Figure |11| 
reveal this isotropy to a high degree, and this is a non¬ 
trivial check of our method, since rotational symmetry 
in the bosonic model is only emergent. It is violated by 
the ultraviolet regularization (see Appendix]^. This fact 
also explains the small differences that are still visible 
between the transverse and longitudinal components. 

In Figure we plot the same data as in Figure [g 
but now in rescaled units, according to (69). We ignored 
small differences between the large distance behavior of 
the transverse and longitudinal components, and scaled 
all components with ^ where the Kondo length 


FIG. 12: (color online) The isotropic correlation functions of 
Figure g in rescaled units 1/^ on the vertical axis and ^ on 
the horizontal axis. The Kondo length ^ was estimated 
from Eq. (1^. 


was calculated using Eq. (57). We clearly see that it is 
possible to scale correlation functions calculated at dif¬ 
ferent J onto universal curves. As expected, we observe 
that the cross-over from slower than 1 /x decay inside the 
cloud to decay at large x occurs around x ^ The 
precise behavior of the Kondo length ^ as a function of 
J will be analyzed further in the next subsection. At 
this point we note that it varies from 11a at J = 0.84 to 
366a at J = 0.42. Such a large variation in the spatial 
scale implies that the observed scaling is non-trivial and 
reflects the universality of the Kondo problem. 


In Figure [T3| we compare the universal scaling functions 
for the Okp and the 2kp components of the cloud, by 
plotting them on top of each other. We determined the 
single universal scaling functions by fitting high order 


polynomials through the scaled data set of Figure 12 


We clearly see that the 2kp component dominates the 
Okp component at small x, consistent with the known 


small X asymptotics of Eq. 70 
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FIG. 13: Universal scaling curves on the isotropic line: the 
black curve represents the longitudinal component Xo, and 
the gray curve represents the transverse component X 2 kp ■ 
These two single curves were obtained by fitting high order 
polynomials through the rescaled data set of FigureThin¬ 
ner straight lines indicate the pure power laws X{x) oc 1/x 
and X{x) oc 1/x^. 


3. Kondo cloud in the anisotropic case 


Having confirmed that the cloud displays the expected 
universal scaling for isotropic couplings, we move on 
to the general anisotropic case. The existence of two 
independent couplings Jy and Jl, or equivalently a 
and A, implies that the universal scaling picture is less 
straightforward. To guide our investigation, let us review 
known results obtained by an improved poor man’s scal¬ 
ing argument, applicable with greater generality in the 
anisotropic case. (Standard poor man’s equations (63) 
can only be trusted when both Kondo exchange couplings 
are small.In the language of the spin-boson model, it 
is known that to leading order in A, but arbitrary a, 
increasing the short-distance scale a by adl is approx¬ 
imately equivalent to changing a to a da and A to 
A + dA, where 


da f A \2 
- = -(oA) a, 


— (aA) = (1 - Q!)(aA). (71) 


Integration of these flow equations yields scaling trajec¬ 
tories 


(«A)^ 


— a + ln(Q;) = constant. 


(72) 


or equivalently, in the language of the Kondo model 


t2 j2 


(27rf 




= constant. 


(73) 


A few of these trajectories in the Jy - J_l plane are plot¬ 
ted in Figure We stress that the standard weak cou¬ 
pling renormalization equations (63) cannot be trusted 


n 



Jl 


FIG. 14: Scaling trajectories that derive from the improved 
poor man’s scaling equations (711, that are correct to leading 
order in J±. In order to emphasize that only the small J± 
part of trajectories are to be trusted, dashed lines are used 
for J±I2'k > 0.1. 



J|| 


FIG. 15: (color online) Map of the parameter space of the 
Kondo model, where dots indicate points where we have col¬ 
lected data to investigate the scaling behavior of the screen¬ 
ing cloud. Below the lower solid curve, Tk ^ 10“®/a and 
we generally find that more than M = 7 coherent states are 
required for a converged result. Above the upper solid curve, 
Tk ^ 1 /u and the screening cloud is poorly resolved because 
of non-universal ultraviolet effects. At the points included in 
the shaded region, we find that the screening cloud is nearly 
isotropic down to distances deep inside the cloud. 


in the regime where Jy ^ 1, and that is why we have to 
work non-perturbatively in a. It is also important to re¬ 
member that these trajectories are only meaningful at J± 
sufficiently smaller than tt. Always bearing this proviso 
in mind, the following statement holds: screening clouds 
associated with two Kondo Hamiltonians, whose parame¬ 
ters lie on the same scaling trajectory, can be scaled onto 
the same nniversal line shape. 

Below, our work will be guided by two qualitative fea¬ 
tures of Figure 14 The first is that the flow trajecto¬ 
ries intersect the Jy-axis at 90°. For large Jy, the parts 
of the trajectories that may reasonably be trusted, run 
nearly parallel to the J± axis. Going beyond the small J_l 
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FIG. 16: Onset of spin-isotropy in the Qkp component of the Kondo cloud, black curves corresponding to the longitudinal 
correlator oIXqI and gray curves to the transverse correlator a|X5^|/2. Top left panel is a = 0.3, top right panel is a = 0.45, 
bottom left panel is a = 0.6, and bottom right panel is a = 0.8. In each panel, different curves of the same shade correspond 
to a few selected values of A. 


regime of the figure, we know that, at the Toulouse point 
J|| = 27r(l — l/-\/2), the scaling trajectory runs exactly 
parallel to the J_l axis. (See Appendix [P}) Universality 
at fixed J|| is also an obvious feature of the exact analyti¬ 
cal expression for the Kondo overlap ( [62| ). For J > 0.457r 
{a < 0.6), we will therefore try to scale clouds at fixed 
J||, but different J_l, onto each other. 

The second pertinent feature of Figure is that for 
^ 7r/2, trajectories are approximately hyperbolas 
J| — Jj_ = c. In other words, the isotropic line Jy = J± 
is an attractor for the renormalization flow. This im¬ 
plies that at sufficiently large distances x > Xiso, the 
cloud must tend to the isotropic cloud. However, for 
X < Xiso, the poor man’s scaling picture of Figure 
indicates pronounced anisotropic behavior. This short 
distance region corresponds to the scales that have to be 
integrated out for an anisotropic point (Jy, J_l) to flow 
close to the isotropic line. Note, furthermore that the 
above discussion only applies to points that are in the 
general vicinity of the isotropic line, and are associated 
with a Kondo temperature sufficiently lower than 1/a. If 
a point is too far from the isotropic line to start with, or 
if the initial Kondo temperature is too large, the size of 
the cloud renormalizes down to the ultraviolet scale, be¬ 


fore the renormalized couplings become isotropic. At this 
point the notion of flow trajectories in a two-dimensional 
parameter space breaks down, and Kondo physics mixes 
with ultraviolet physics. Based on these observations, we 
will compare the cloud calculated at Jy, J_l ^ tt/S (i.e. 
reasonably close to the isotropic line and with a decent 
Kondo length) with the isotropic cloud. 

We start the presentation of our results in the regime of 
small a (large Jy). This regime does not overlap with the 
isotropic Kondo regime, because the size of the screen¬ 
ing cloud flows to the short distance cut-off a before the 
renormalized couplings come close to the isotropic line. 
For small a, the single coherent state (Silbey-Harris) 
approximation is accurate (see Figure J3, and we have 
derived simple analytical expressions ( [5^ for the cloud 
correlation functions, that reveal the following universal¬ 
ity. By appropriately rescaling the four cloud compo¬ 
nents and the coordinate x, clouds calculated at different 
J_L <C TT but fixed Jy (i.e. fixed a) have the same line 

shape: For the longitudinal components Xj , j G {0, 2kp} 
one defines 

Xjia,x) = ^llXj{x), i = a;/Cy, (74) 

and finds from the Silbey-Harris result the simple ana- 
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and transverse cloud components plotted against the scaled 
distance Perfect isotropy corresponds to 2Xq/Xq — 1. 

Solid curves correspond to a = 0.85 and various values of A 
ranging from 0.156/o to 0.305/a. Dashed curves correspond 
to Q = 0.75 and various values of A ranging from 0.120/a to 
0.268/a. The Kondo length was calculated for each curve 
from Eq. (|57|). For the set of data at a = 0.85, varies from 
18.0a to l58a, while for a = 0.75, varies from 12.6a to 
108a. 


lytical formulas: 


Xq (a, x) = —ReF 




yax\ 


- )’ 


2y/almF 

^i^/ax\ 

V - )\ 


(75) 


In the small a regime, the two transverse components 
Xj-, j e {0, 2fcF}, are equal. One defines 


Xfia,x) = ilXj{x)/U, X = x/^||, 


(76) 


and finds the analytical expression: 


XHa,x) = - —exp 


—2y^Iie F 


i^OLX 


(77) 


The expression for Xg (x) could suggest a further rescal¬ 
ing X = ^Jax that would lead to an a independent line 
shape for Xg. However, one cannot simultaneously get 
rid of the a dependence in the other components of the 
cloud. Furthermore, it turns out that X^\a,z/^/Q) is no 
longer a-independent beyond a > 0.4 (as shown by our 
numerical simulations below). 

The two pertinent facts to emerge from this discussion 
are the following. Firstly, whereas the magnitudes of 
the longitudinal and transverse components of the cloud 
have to be rescaled by different inverse lengths l/^y and 
C-l /(?||)^5 tlie coordinate x is always rescaled by In 
other words, the longitudinal and transverse clouds have 
different characteristic magnitudes l/^y and ^_l/(^||)^, 


but the same characteristic size ^y. Secondly, the uni¬ 
versal scaling curves depend in general on a. Thus, as 
expected from poor man’s scaling, clouds at the same Jy 
but different Jl have the same shape, whereas clouds at 
different Jy in general have different shapes. 

The differences in scaling behavior between the small 
a and isotropic clouds are thrown into sharp relief when 
one considers the small x asymptotic form of the cloud 
correlators. At a; <C 1, the small a results (75) and ([77]) 
lead to asymptotic behavior: 


X^ = X. 


2kf 




Xg ~ Ini, X. 


‘2k F 


— X 


(78) 


The Xg component diverges much more slowly in the 
i —>■ 0 limit than is the case at isotropic couplings 
(Ini vs. —l/i(lni)^), whereas the x|^^ component di¬ 
verges more rapidly than in the isotropic case (—1/i vs. 
1/ilni). Thus, the alternation between ferro- and an¬ 
tiferromagnetic correlations close to the impurity is en¬ 
hanced in the longitudinal direction. This is not surpris¬ 
ing, as small a implies Jy 3> J_l. Furthermore we see 
that the small x asymptotic behavior of the transverse 
cloud is explicitly a-dependent, displaying power-law be¬ 
havior with an exponent —2(1 — y/a) = — Jy/yr, implying 
a divergence as i —>■ 0 that is closer to 1/x^, than to 1/x. 

The next obvious question is how the small a scaling 
picture evolves as a increases. In order to investigate 
this, we calculated the cloud correlators at the points in 
parameter space shown in Figure |15[ The most obvi¬ 
ous feature of the data, is that the cloud becomes more 
and m ore isotropic as a increases. This is revealed in 
Figure 16 where we plot xj(a;) and X,^(a;)/2 on top of 
each other for respectively a = 0.3, a = 0.45, a = 0.6, 
and a = 0.8, in separate panels. Each panel contains 
curves for several values of Jj_. At a = 0.3, the trans¬ 
verse and longitudinal clouds are very different. As a 
increases, the differences become smaller. By the time 
we reach a = 0.8, the cloud is isotropic to a high degree 
of accuracy for all values of J_l considered. The behav¬ 
ior of the x|j,^ (x) and X^^ (x) components of the cloud 
(not shown) present very similar behavior. We find simi¬ 
larly isotropic clouds at all the points in Figure [Ts] inside 
the shaded region. 

Poor man’s scaling suggests that for anisotropic ex¬ 
change couplings, there is a scale Xiso below which 
isotropy breaks down. To investigate this, we plot the 
ratio 2 Xq/X 5^ as a function of a:/^y (Figure 


17). 


We do 

so for the data collected at a = 0.75 and for cT^ 0.85. In 
both cases we see that 2 Xq / X^ reaches a plateau that is 
within 10% of unity. This is the same degree of isotropy 
as we obtain on the isotropic line itself. In both cases, 
the plateau is reached at a distance that is a small frac¬ 
tion of the Kondo length ^_l. In other words, isotropy 
extends deep inside the cloud. The ratio Xiso/Cy between 
the scale where isotropy sets in, and the Kondo length, is 
more or less constant at fixed a, while the Kondo length 
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FIG. 18: Scaling curves for the four component of the Kondo cloud computed at hxed a = 0.6 (Jy = 1.42). Raw data is shown 
in the insets, with different curves correspond to eight different values of A collected in the interval A € [0.05/a, 0.25/a] (i.e. 
Jj_ € [0.157,0.785]). Excellent scaling behavior is obtained for all components. 


itself varies by an order of magnitude. From the figure 
one can read off that Xiso/C|| ~ 0-2 for a = 0.75 and 
2^iso/C|| 0-05 for a = 0.85. 

Next, we investigate the precise scaling behavior of the 
data. For a < 0.6, we take our cue from the small a re¬ 
sults we obtained above, together with the understand¬ 
ing from poor man’s scaling, and we thus scale corre¬ 
lation functions according to the Ansatz (74) and (76). 


In Fig. we show raw and scaled results for a = 0.6. 
Similar figures for a = 0.3 and 0.45 can be found in 
the supplementary material.l^ For all values of a that 
we considered, the data nicely scale onto single universal 
curves. We conclude that for practical purposes, scaling 
Ansatz of the form (74) and (76) hold for all values of 
a < 0 . 6 . 

These universal scaling curves are a-dependent and dif¬ 
fer from the isotropic scaling curves. This can clearly be 
seen in Figure |19[ where we compare the universal line 
shapes of the scaled cloud at a = 0.3 (top panel) and 
at a = 0.6 (bottom panel) to the universal isotropic line 
shapes obtained in Sec. |VIII D 2[ Although the single 
coherent state (Silbey-Harris) approximation underesti¬ 
mates the Kondo length by 50% at a = 0.3, the simple 
equations ( 75||77 ) still predict the universal line shape of 
the screening cloud very well. Inside the cloud (a: < ^y). 


these universal curves are very different from the univer¬ 
sal isotropic curves. As a increases, deviations from the 
Silbey-Harris line-shape set in, and scaling curves start 
to resemble the isotropic curves more closely, as can be 
seen for a = 0.6 in the bottom panel of Figure [T^ 

A prominent feature of the anisotropic scaling curves 
is that, at small distances, there is a stronger suppression 
of the Ag component over the X^f.^ component, than in 
the isotropic case. As a result, the alternation from ferro- 
to antiferromagnetic correlations between the impurity 
and the electron gas inside the cloud is enhanced in the 
longitudinal direction. This is illustrated in Figure [20} 
where we compare the full longitudinal cloud X^^{x) = 
Ag(a:) -I- cos{2kFx)X2j.^{x) calculated at a = 0.3 with 
the result on the isotropic line. 

To complete the scaling analysis, we consider the 
screening cloud at 0.75 < a < 0.9 (the shaded region in 
Figure [T^. In this region, we have seen that the cloud is 
approximately isotropic. In Figure |21[ we therefore plot 
the transverse and longitudinal components of the cloud 
on the same graph. Insets show raw results in units of a 
on the horizontal and 1/a on the vertical axes. The main 
panels show data in rescaled units of ^y on the horizontal 
and 1/^y on the vertical axes. We see that the scaled 
data collapse very well onto single scaling curves. As be- 
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FIG. 19: (color online) Comparison between the universal 
line shapes of the scaled cloud (symbols in color), at q = 0.3 
(top panel) and at a = 0.6 (bottom panel), to the universal 
isotropic line shapes (full and dashed black lines). The curves 
were obtained by fitting data sets such as those displayed in 
the main panels of Fig.[^to polynomials of high degree. On 
the vertical axis, the longitudinal components are plotted in 
units of 1 /5||, while the transverse components are plotted in 
units of ■ The isotropic curves are those of Figure 
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fore, the Kondo length was calculated using Eq. (57), 
i.e. no fitting parameters were used to obtain the high 
degree of collapse. It is instructive to compare the main 
panels of Figure for all the data collected, to Figure 
[T^ that shows scaled data only for the subset of points 
that lie on the isotropic line J_l = J\\. We see the same 
degree of collapse onto single curves in both figures. In 
other words, within the numerical accuracy of our calcu¬ 
lation, the universal line shape of the screening cloud is 
the same for all points in the shaded region of Figure [TS] 
To further confirm this conclusion, we fitted a high order 
polynomial through the complete scaled data sets of Fig¬ 
ure In Figure [2^ we compare this ht to the universal 
scaling curves we obtained previously by only considering 
the cloud on the isotropic line (Figure [I^. We see nearly 
perfect agreement. Of course, the screening cloud at an 
anisotropic value of the exchange couplings, Jj_ ^ Jy, 
only follows this universal line shape for distances larger 





FIG. 20: The full longitudinal cloud X''\x) = 7fg(a;) -|- 
cos{2kFx)xk (x) at a = 0.3 (black) and on the isotropic 
line (gray), as shown by full lines (dashed lines are the re¬ 
spective envelope functions 7fg(2;) ± (x)). Both curves 

were normalized using their magnitude at a; = 2^Ikp. For the 
purpose of showing the oscillations, we took kp = 207r/^x. 


than Xiso, the scale at which isotropy sets in. However, 
as we have seen in Figure there is a sizable region 
around the isotropic line where isotropy already sets in 
deep inside the cloud. 

Finally, we plot in Figure 23 the Kondo length ^y, 
which corresponds to the characteristic size of the cloud, 
for some of the data we have collected. The large range 
over which ^y varies for each value of a, confirms that the 
scaling behavior we see is non-trivial. We compare the 
calculated value of ^y to ^ oc 1/Tk, with Tk the standard 
poor man’s scaling estimate (64) of the Kondo tempera¬ 
ture. With a logarithmic scale on the vertical axis, which 
hides mismatches by factors of order unity, we find good 
agreement for the isotropic and the a = 0.8 data sets. 
That these sets yield the best agreement is expected, as 
they are closer to the point Jx = Jy = 0 in parameter 
space than the other sets, and the version of poor man’s 
scaling that leads to the estimate (64) for Tk assumes 
small exchange couplings. 


IX. SUMMARY OF RESULTS AND 
CONCLUSION 

This paper has provided a very extensive study of the 
Kondo screening cloud for a wide range of parameter 
values, including the spin-isotropic and strongly spin- 
anisotropic regimes. Methodologically, we derived simple 
but controlled analytical expressions in the case of strong 
anisotropy, and developed an original and very powerful 
numerical technique to tackle the problem in its com¬ 
plete generality. This allowed us to investigate the Okp 
and 2kF component of the Kondo cloud correlator sepa¬ 
rately, both for the longitudinal and transverse response. 
In addition, we have examined, again both analytically 
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FIG. 21: (color online) The screening cloud calculated at all 
points inside the shaded region (0.75 < a < 0.9) in Figure [T5| 
Insets show the raw data, while the main panels show data 
in units of on the horizontal axis and units of l/^y on the 
vertical axis. The top panel contains the transverse and longi¬ 
tudinal 0 /cf components plotted together. The bottom panel 
contains the transverse and longitudinal 2kp components. 


and numerically, the structure of Kondo overlaps intro¬ 
duced in recent works 

Our main results concern the universal scaling of corre¬ 
lations between the impurity spin and the electron spin 
density in the anisotropic Kondo model. They can be 
summarized as follows. At large a;, transverse correlators 
equal X^{x) = while longitudinal correlators 

equal x\{x) = —(here k = Q, 2kp refer to respec¬ 
tively the forward or backscattering component of the 
correlator). In general, the two emergent length scales 
^_L and ^11 are not equal, and we find the following uni¬ 
versal scaling behavior: 

Xt{x) = ^Xi{a,x/(,\\), 

X|(x) = llf(a,x/^l|). (79) 

The scaling curves are Jj^-independent, but remain a- 
dependent (and therefore Jy-dependent). As can be seen 



FIG. 22: The solid gray curves represent the (isotropic) uni¬ 
versal scaling curves inferred using all the data of Figure 
i.e. clouds calculated at all the points inside the shaded region 
in Figure [T^ These curves were obtained as the best fit of a 
high order polinomial through the scaled data of Figure |21| 
The dashed black curves are the universal scaling curves in¬ 
ferred using only clouds calculated at perfectly isotropic cou¬ 
plings (See Figure [l3|). 



J± 

FIG. 23: (color online) The Kondo length as a function 
of J± for several a values. The solid lines shows ^ = c/Tk, 
with Tk the standard poor man’s scaling estimate ( |64| | and 
c a fitting parameter. The same value of c = 0.3 we used for 
all five curves. 


from the scaling equations, the two lengths ^_l and ^y 
have fundamentally different roles. The length ^y sets the 
size of the screening cloud, whereas the ratio C-l/C|| sets 
the relative magnitude of transverse correlations, com¬ 
pared to longitudinal correlations. 

As Jy decreases, the explicit Jy dependence of the scal¬ 
ing curves becomes weaker and weaker. For Jy ^ 1, 
we find no discernible Jy -dependence anymore, and the 
universal scaling curves become those of the isotropic 
model. In this region of parameter space, we find that 
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isotropy already sets in deep inside the cloud. We man¬ 
aged to calculate the universal scaling curves down to 
distances = 0.1 or less. This far inside the cloud, 
the backscattering {2kp) component of the longitudinal 
correlator dominates the forward scattering (Okp) com¬ 
ponent. Inside the cloud, correlations between the impu¬ 
rity spin and the electron spin density therefore alternate 
with wave vector 2kp between being ferro- and antifer¬ 
romagnetic. The effect is more pronounced in the regime 
J|| S> J_L, but we can still clearly resolve it when the 
cloud becomes isotropic. 

In previous numerical studies, results were obtained 
using the numerical renormalization group (NRG). This 
required solving a two-impurity problem, where the po¬ 
sition of the second (fictitious) impurity represents the 
point X where the correlator is evaluated. An indepen¬ 
dent two channel Kondo problem must be solved by NRG 
at every position where the correlator is evaluated, mak¬ 
ing the calculation computationally expensive. In addi¬ 
tion, it proved difficult to resolve in NRG the dominance 
of the backscattering components of the cloud over the 
forward scattering components inside the cloud. 

We have shown that the coherent state expansion, the 
method we employed in this Article, complements pre¬ 
vious numerical renormalization group studies, and even 
offers simple analytical insights in the regime of strong 
spin-anisotropy. Thanks to a technical improvement that 
dramatically reduces the number of variational parame¬ 
ters, it offers excellent accuracy at a reduced computa¬ 
tional cost, and does not require one to discretize the 
bath degrees of freedom. Yet, we note that the numeri¬ 
cal renormalization group was successfully used to study 
finite temperatur^ and time-dependeniP features of the 
Kondo cloud. We envisage future work to extend the 
coherent state expansion into these domains as well. 


Appendix A: Mapping between the Kondo and 
Spin-Boson Models 


We separate spin and charge degrees of freedom by 
defining: 

ag = {bg^ + bgi)(V2, bg = {bg^ - bgi)jV2. (A3) 

We also define new SU(2) operators, in which fermion 
number and impurity degrees of freedom are mixed. 

g- = F^Fia~, s+ = (s“)'l', Sz = cTz- (A4) 

In the bosonic representation, the Hamiltonian reads 


Ho ^J2q(alag + blbg), 


q>0 


/I'v.n V 


q>0 

H 

^ 2'Ka 1 




(AS) 


where 


(^t(0)-</>t(0) 

V2 


= (A6) 

q>0 * ^ 


We note that spin {bg) and charge {ug) degrees of free¬ 
dom are decoupled, and that only spin degrees of freedom 
couple to the impurity. We restrict the Hamiltonian to 
the vacuum of the charge sector and drop the a^a terms 
in future expressions. 

Finally, we bring the Hamiltonian into the standard 
form of the (unbiased) spin-boson model, via a unitary 
transform 


U = eicp{—iSz4’/V2). (A7) 


Under the action of U, 


UbgU^ = bg 



Us-u'f = 

(AS) 


In this Appendix, we review the mapping from the 
Kondo Hamiltonian @ to the spin-boson Hamiltonian 
(§. The first step is to express the (even) fermionic de¬ 
grees of freedom in terms of bosonic ones bga-, q = 2Tm/L^ 
n G {1, 2, 3, ...} with [bgo,b^g,„,] = 5gg>5aa'- For this 
purpose, we invoke the bosonization identities 


and 


Phr 

bq<y = \ 1- dxe~^^^il}l{x)i)a{x), 

V Lq J_p /2 


Fl 


(Al) 


,l;Ux) = 


,{x) = + h.c.. 

q>0 * ^ 


(A2) 


Here F^ denotes the unitary Klein factor operator asso¬ 
ciated with increasing the number of a electrons by one. 


so that UHW = Hsb + const, where Hsb is given by ([^. 


Appendix B: Correlation functions in the bosonic 
representation 


In this Appendix, we transform the fermionic cloud 
correlators 0 onto their bosonic counterparts ( [To| ). To 
do so, we represent the electron spin density in terms 
of bosonic degrees of freedom, a nd s ubsequently perform 
the unitary transformation of (A7) on these operators. 
To this end, it is useful to note from 0 the following 
identities: 


1 

71 

1 

71 


i>a{x) F-lpai-x) 
-Ipaix) - i^ai-x) 


(Bl) 
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The ‘approximately equal’ signs in these equations in¬ 
dicate that the operators on the right hand side repre¬ 
sent versions of the point operators on the left, that have 
been broadened by an amount ^ a, inversely proportional 
to the ultraviolet cut-off. Now consider the expectation 
value which, from Q, can also be written as 

X-^(x) = 4Re ^cr"'0|(a;)Vi^(x)^^ . (B2) 


Because the ground state has good spatial parity, X-*- (x) 
is an even function, and we may write 


X-‘~(x) = 2Re ^cr 




The last expression is manipulated to become 



X-^(x) = Re < (cr iphx) + tpl(-x) ip^(x) + ip^(-x) 


J / K 


+ (cr 


ip\{x) - ipli-x)^ - ipii-x)^ } . (B4) 


=0 


The term marked with the underbrace is zero because of the following reason: 'i/'|(x) — x) creates an electron in 
an odd parity single particle orbital, and these orbitals do not couple to the impurity, so that: 

(|cr" ^/;|(x)-^^|(-x) ^jji{x)-^i{-x) = ^{x)-'tPl{-x) '^i{x)-^i{-x) (B5) 


Sinc e the antiferromagnetic Kondo ground state is a singlet, (cr )j^ = 0. Substituting from ( |Bl[ ) into the first line of 
(B41, we obtain 

X-*-(x) = 2Re (^cr“ tp\{x)ilJi{x)+tj;^-x)ilJi{-x) + e~‘^'‘’'‘^''ilj}.{x)'tpi{-x) + e^''''’"^tp}-{-x)tpi{x) . (B6) 

Performing a similar calculation for Xll(x), we find 

[ipl{x)^a{x) + 'iIjI{-x)jP^{-x) + e-^^'^’='^^l{x)'il;^{-x) + (B7) 


where sgn(t) = -I- and sgn(|) = —. We finally manipu¬ 
late the above expressions as follows: 


1. The fermion operators are bosonized. 


2. We apply the unitary transformation that maps 
the Kondo Hamiltonian onto the spin-boson model. 
In the resulting expressions, expectation values are 
with respect to the ground state of the spin-boson 
model, indicated without a subscript. 


3. We normal order the bosonic operators. In the re¬ 
sulting expressions, the operators oj and a,, are 
replaced by zeros, since annihilates the ground 
state. 


Finally, we take the large system limit, which allows us 
to make replacements such as I — 27r(a-|- 

ix)jL. The resulting expressions are those presented in 
the main text (101. 


Appendix C: Contour integral for cloud calculation 


The integrals involved in calculating the cloud are of 
the form: 


dqe 


-(a±ix)q 


(9), 


(Cl) 


where the function r{q) goes to zero as least as fast as 
1/q for large q. Furthermore, the only non-analyticities 
that r[q) possesses in the complex q plane are (first or¬ 
der) poles at Wn, n = 1, ..., M. These integrals can be 


performed by extending the method employed in Sec. VI 
where the energy was calculated. 

Define z = a±ix and 9 = —arg(z). Integrating around 
the shaded area A, whose azimuthal boundary is under¬ 
stood to lie at |g| —00 in Figure 24 we obtain 


/ dge-^«r(9) = e*® / dpe-I^IPr(pe*®) 

/o Jo 


2TTi ^ Res 


"(9)d 


(C2) 


uj„eA 


The integral on the right hand side on the first line is 


of precisely the same form as those considered in Sec. VI 
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We also note that the summation on the second line van¬ 
ishes if there are no poles ujn with positive real parts. 
Using the results of Sec. VI we obtain: 


N 

dqe~^‘^r{q) = ^ Res [e~^‘‘Ei{-zq)r{q),uJn] 

m—1 

+ 2-Ki ^ Res [e~^‘^r{q),ujn] . (C3) 



Appendix D: Exact results at the Toulouse point 


In this appendix, we review the mapping between 
the Kondo model at the Toulouse point a = 1/2 (i.e. 
J|| = 27r(l — l/-\/2)), and the non-interacting fermionic 
resonant level model, following Ref. [43l We also derive 
an exact expression for the longitudinal Okp component 
of the screening cloud. The derivation takes the same 
route as that set out in Refs. H^andliHl 


The starting point is the bosonic representation (A5) 


of the Kondo hamiltonian. The fermionic non-interacting 
resonant level i?rim is obtained, by applying the unitary 
transform 7?rim = U'H{U')\ where 


U' = exp 


2 V V 2 J V2_ 


(Dl) 


The bosonic field (j) is defined in Eq. ( |A6[ ), and the oper¬ 
ator 


N_ = {N^ - iVj/2, 


(D2) 


is the z-component of the total spin of the electron gas. 
New degrees of freedom 


= exp 


2Tri f 1 . . , / N 





c-ikx -a\k\/2-^ 


= 


exp 


ZTT I V_ - 


(D3) 


emerge, that obey the usual fermionic commutation re¬ 
lations. Here 


(/(x) = 


4)^{x) - ^;(x) 


(D4) 


with defined in Eq. (A2|, and denote Klein 

factors. In terms of these new fermions, i?rim reads 


^Irim = qclcq + d + d^^{x)] . (D5) 

^ ' V 27ra 


At the Toulouse point, the unitary transformation U' 
maps the ground state of the Kondo Hamiltonian onto 
the ground state of i?rim- Since Hrim is a quadratic 
Hamiltonian, ground state correlation functions can be 
calculated straight-forwardly for transformed observ¬ 
ables, provided they are expressed simply in terms of the 
new fermionic degrees of freedom. Here we focus on the 
longitudinal Qkp component Xq for which this is easily 
accomplished. 

Applying the unitary mapping to the z-components of 
the impurity spin and the electron spin density operators, 
one finds 


^0 (^) =v^ {^^dU- 0 [i/j(x)'>'^(x) + ip(-x)'<'i/;(-x)] ^ 
-(2-V2)^^, (D6) 

where the expectation value on the right hand side of the 
first line is with respect to the ground state of i?rim- This 
expression is further simplified using Wick’s theorem and 
noting that the occupation probability for the resonant 
level is one half, i.e. = 1/2. In terms of the orig¬ 

inal degrees of freedom, this corresponds to (az) = 0 
which is guaranteed by the singlet nature of the Kondo 
ground state. Working directly with idrim, this follows 
from particle-hole symmetry. (Note that the energy of 
the resonant level co-incides with the Fermi energy.) To 
be more precise, Idrim is invariant under the combined 
action of particle-hole and spatial inversion. Owing to 
the same symmetry, (dt'(/(x)) = (dt'i/)—x)) . Putting 
everything together, one obtains 

Xll(x) = -2^2 |(dt^(x))f - (2 - V2)^^. (D7) 

The expectation value (dt^/(x)) can be calculated from 
the known single particle Green’s functions of the reso¬ 
nant level model. In the thermodynamic limit (L —>■ 00 ) 
the exact result is 

(dii/Cx)) 


J± 


di?(x) 

dR{-x) 

y/2'Ka j 

'-00 Stt 

Lo — S(a;) 

oj - E(tLi)*_ 
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where 0r{x) is a regularized step function 

C°° ^i{k-ui)x-a\k\ 


Or{x) =i f 
J —( 


0{x) + 


27r a;+ — k 

jr(0, (jj{ix + a)) 


27r 


„-u,a *r(0,w(ia:-a)) 


27r 


(D9) 


with w"*" = w + iO"*", and S(a;) is the retarded resonant 
level self-energy 


j? 




g-a|fc| 

W+ - fc’ 


[e“‘^r(0, aw+) - e-“‘^r(0, -aw+ll . (DIO) 

dTT^a 


In general, the integral (D8| has to be evaluated nu¬ 


merically. However, for |a;| ^ a, the following approxi¬ 
mations are valid, which yields an analytical result. As 
a function of (negative) w, the regularized step function 
9r{x) is su ppre ssed as e‘^“ for a; > 0. As a result, the in¬ 
tegrand of (D8) is strongly suppressed for — w > 1/a. As 
a function of x, the regularized step function is smoothed 
at a scale \x\ < a. For —w < 1/a and |x| S> a, the regu¬ 
larized step function can be replaced with the sharp step 
function 0[x). For x ^ a, the rapidly oscillating factor 
^luix qq: integral at w ~ 1/a; <C 1/x, so that 

the self-energy S(a;), that varies on the scale of 1/a, can 
be evaluated at w = 0 


7 P 
47ra 


(Dll) 


For |x| 3> a, one then finds 

J± 


{d)ip{x)) = 


dio 


j\x\ 


y/PPci J —oo m ■ 


■ 

^ 47ra 


J± 


( 

I 47ra J 


F - 


\/27ra 


2tt 


(D12) 


where F{z) is defined in Eq. (38). We substitute this 


into ( |D7[ ) and also drop the second term in that equation, 
which is valid for |x| a, to obtain 


^ TT^ 47ra V 47ra 


(D13) 


The Kondo length is defined as in (56) in the main 
text, i.e. ^11 = lim 2 ;_>.oo x^Xg (x). This yields 


^ ^/2 Aira 

^11 


and a universal scaling function (cf. Eq. (@) 


x!{a = l/2,i) = -^Fl-^"j . (D15) 
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Scaling analysis of the anisotropic cloud 


In the main text, we showed extensive data at a = 0.6 (i.e. Jy = 1.42) demonstrating that cloud correlation 
functions calculated at different values of A = Jx/ttc scaled onto the same universal functions (Figure 18). Here, 
we show for completeness data that demonstrate scaling at a = 0.3 (Figure [sTj), a = 0.45 (Figure [S^, and a = 0.8 
(Figure [S|). 






FIG. SI: Scaling curves for the four components of the Kondo cloud computed at fixed a = 0.3 (Jy = 2.84). Raw data is shown 
in the insets, with different curves correspond to five different values of A collected in the interval A G [0.005/a, 0.07/a] (i.e. 
J± G [0.0157,0.220]). 
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FIG. S2: Scaling curves for the four components of the Kondo cloud computed at fixed a = 0.45 (Jy = 2.07). Raw data is 
shown in the insets, with different curves correspond to six different values of A collected in the interval A G [0.01/a, 0.125/a] 
(i.e. Ju G [0.0314,0.393]). 
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FIG. S3: Scaling curves for the four component of the Kondo cloud computed at fixed a = 0.8 {J\\ — 0.663). Raw data is 
shown in the insets, with different curves correspond to eight different values of A collected in the interval A G [0.12/a, 0.28/aj 
(i.e. J± G [0.377,0.880]). 


























